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Abstract. A novel tool for tsunami wave modelling is presented. This tool has the 
potential of being used for operational purposes: indeed, the numerical code VOLNA is 
able to handle the complete life-cycle of a tsunami (generation, propagation and run-up 
along the coast). The algorithm works on unstructured triangular meshes and thus can 
be run in arbitrary complex domains. This paper contains the detailed description of the 
finite volume scheme implemented in the code. The numerical treatment of the wet/dry 
transition is explained. This point is crucial for accurate run- up/run-down computations. 
Most existing tsunami codes use semi-empirical techniques at this stage, which are not al- 
ways sufficient for tsunami hazard mitigation. Indeed the decision to evacuate inhabitants 
is based on inundation maps which are produced with this type of numerical tools. We 
present several realistic test cases that partially validate our algorithm. Comparisons with 
analytical solutions and experimental data are performed. Finally the main conclusions 
are outlined and the perspectives for future research presented. 
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1. Introduction 



After the 2004 Boxing Day tsunami [SB06] and the 2011 Honshu, Japan tsunami, there 
is no need to explain the importance of research on tsunami waves. One of the primary 
objectives in this field consists in establishing and developing Tsunami Warning Systems 
(TWS) [Tat97, TGB+05] and inundation maps. This task is non trivial as explained by 
Synolakis [Syn05]: 

For reference, the United States and Japan took more than 20 years to 
develop validated numerical models to predict tsunami evolution. And it 
took the US National Oceanic and Atmospheric Administration 30 years 
to fully develop its bottom-pressure recorders, which have been reliably 
detecting tsunamis for the past ten years. 

After the Boxing Day tsunami, while developing their own national and regional capabili- 
ties, countries in the Indian Ocean and the Caribbean Sea have asked the PTWC (Pacific 
Tsunami Warning Center) to act as their interim warning center. India and Australia 
now have fully working national centers, while the National Oceanic and Atmospheric 
Administration of the U.S. has assisted both with instrumentation and the sophisticated 
forecast technology used in the Pacific. Europe however is trying to reinvent the early 
warning wheel. As a result, the Mediterranean remains the only world sea unprotected by 
any warning system. The 2011 Honshu tsunami also showed that the tsunami community 
did not do enough to anticipate future events, even in Japan which is arguably the most 
tsunami-ready nation in the world. 

The mathematical modelling and computation of propagating tsunami waves play an 
important role in TWS. Precision and robustness of the algorithm will affect performance 
and reliability of the whole system. 

The importance of tsunami generation modelling is often underestimated by the scientific 
community. During several years the research of our group was focused on this topic 
and interesting results were obtained [DD07c, Dut07, DDK06, KDD07, DD09b, DDIO, 
DMGDIO, DMCSIO]. We tried to incorporate some recent developments [DD07c] from 
this field into the VOLNA code. 

The recent events in Japan should convince the scientific community of the urgency 
to complete inundation maps. These are maps that show the extent of possible tsunami 
flooding from hypothetical earthquakes in their vicinity. Even in the U.S., only California 
has completed its mapping efforts. Alaska and Hawaii, the most vulnerable U.S. states, do 
not have modern tsunami flood maps for all their coastal communities.^ 

""^Whcn will we learn, Newsweek, March 13, 2011, by Costas Synolakis 
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It is difficult to find a topic in numerical analysis of hyperbolic PDEs which has been stud- 
ied more than the numerical solution to the Nonlinear Shallow Water Equations (NSWE). 
The numerical scheme presented in this paper is not completely novel. The discretization 
methods used in VOLNA can be found in the modern literature on finite volumes methods 
[Kro97, BO04]. The main purpose here is to present a tool for tsunami wave modelling 
which covers the whole spectrum from generation to inundation. The emphasis is on the 
technical work which is typical of a numerical analyst and software developer. Tsunami 
practicioners can then concentrate on the physical aspects of tsunami propagation. 

Nowadays, one is facing a somewhat strange situation. On one hand, there are only a 
few truly operational codes for tsunami wave modelling: MOST, NAMI, ComCot [Ima96, 
TG97, GOSI97, LWC98]. The numerical schemes used in these codes essentially correspond 
to the state of the art of the eighties. On the other hand, there is a plethora of NSWE 
codes developed in academic environments [Gla88, Cas90, Tor92, BDDV98, AC99, VC99, 
AGN05, BQ06, GNVCOO, GHS03, ABB+04, KCY07, Geo06, GL06, ZCIM02, CIM+00, 
NPPN06, CFGR+05, WMC06, Geo08, DKK08, CBB06, CBB07] - see [Bar04] for a review 
of nonlinear shallow water theories for coastal waves. These codes use modern numerical 
methods but most of them have not been developed to satisfy the needs of tsunami oper- 
ational research. This is why we had the idea to develop VOLNA . We tried to combine 
modern numerical techniques for hyperbolic systems with real world application-oriented 
design. The VOLNA code can be run efficiently in realistic environments. It was shown 
that natural coasts tend to have fractal forms [SBG04]. Hence, unstructured meshes are 
a natural choice in this type of situations. The ffist demonstration of the applicability 
of VOLNA to real situations was provided by Poncet et al. [PCD+IO] in their study of 
tsunamis generated by landslides in the St. Lawrence estuary. 

The paper is organised as follows. In Section 2 the physical context of the study and the 
motivation for the choice of the mathematical model are presented. Section 3 contains a 
detailed description of the numerical method implemented in the VOLNA code. In Section 
4 we show computations which validate and illustrate the capabilities of VOLNA . Finally 
the main conclusions are outlined and the perspectives for future research presented in 
Section 5. 



In this study we focus on long wave propagation over realistic bathymetry. A sketch of 
the physical problem under investigation is given on Figure 1. Let us explain the main 
assumptions and the domain of applicability of the VOLNA code. 

First we introduce some characteristic lengths. We denote by the typical wave am- 
plitude, by Hq the average depth and by £ the characteristic wave length. Several dimen- 
sionless numbers can be built from these three quantities, but traditionally one introduces 
the following two: 



2. Physical context and mathematical model 
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Figure 1 . Sketch of the fluid domain. 



The first parameter e measures the wave nonhnearity {e <^ 1 means than nonhnearity is 
weak) while the second parameter /x^ quantifies the importance of dispersive effects (yU ^ 1 
means than dispersion is weak). Taking a typical megatsunami offshore with roughly 

ao — 0.5 m, /^o — 4 km, t ~ 100 km 

yields e = 1.25 x 10~^ for the nonhnearity parameter and /x^ = 1.6 x 10~^ for the dispersion 
parameter. Both are weak. Using asymptotic expansions in the small parameters e <^ 1 and 
/i^ ^ 1, one can derive Serre-type equations [Ser53, Per67, MBS03, Lyn06, DD07a, LB09, 
DM10, CLMIO, BCL+11, CDll, CCll, DCMMll]. The effect of dispersion on tsunamis 
has been investigated recently [DT07, IAK"''07, MFS08]. Due to frequency dispersion, 
longer and higher waves travel faster and separate from the shorter and smaller waves, 
leading to a decrease of tsunami height. It is close to the shore that dispersion might play 
a role. Here the short waves may have a local additional effect on wave impact on coastal 
structures, but they hardly play any role for the runup and inundation caused by the main 
and much longer tsunami. Therefore the consequences of neglecting dispersive effects are 
probably not very important from a practical point of view. Moreover no operational codes 
based on the Serre equations exist as of today. The GEOWAVE code, which combines 
TOPICS and FUNWAVE, is still too expensive from a computational point of view to 
be used as a truly operational code. Neglecting the dispersive effects yields the classical 
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Nonlinear Shallow Water Equations (NSWE): 

Ht + V ■{Hu)=0, (2.1) 

{Hu)t + V ■ (^Hu O M + ^H^^ = gHVh, (2.2) 

where H = h + r] is the total water depth and u = {u,v){x,t) is the depth-averaged 
horizontal velocity. Traditionally, g denotes the acceleration due to the gravity and h{x,t) 
describes the bathymetry. 

Remark 1. The bathymetry h{x,t) is allowed to be time- dependent. It is important for 
the problem of tsunami generation by underwater earthquakes, submarine landslides, etc. 
The coupling with seismology is done through this function. Namely, various simplified 
earthquake models [DD07c, KDD07, DD09b, DDIO, DMGDIO, DMCSIO] provide the seabed 
displacements which are then transmitted to the ocean layer. 

In this study, the NSWE (2.1) and (2.2) are chosen to model tsunami generation, propa- 
gation and run-up/run-down. It is computationally advantageous to have a uniform model 
for all stages of tsunami life since many technical problems are thus avoided. The validity 
of the NSWE for tsunami generation was already examined in our previous study [KDD07], 
where an excellent performance of this model was shown for nondispersive long waves. In 
the present paper we show in Section 4 the ability of the NSWE to model the run-up/run- 
down process. For this purpose, comparisons with a laboratory experiment are performed. 
Thus, the chosen complete approach to tsunami wave modelling is very attractive from 
both the operational and research viewpoints. 

The governing equations (2.1) and (2.2) have nice mathematical properties. In partic- 
ular, this system is strictly hyperbolic provided that H > 0. This property will be used 
extensively in the construction of the numerical scheme (see Section 3). 

Let us discuss the eigensystem of the advective flux. First, we introduce conservative 
variables and rewrite the governing equations as a system of conservation laws: 

^ + V-J'iw) = Siw), (2.3) 

where the following notation was introduced: 

w{x, t) : X M+ ^ w = {wi, W2, Ws) = (H, Hu, Hv), 

/ W2 \ 



Hu Hv 
T{w) = \Hu'^ + f iJ2 Huv 

Huv Hv^ + ^H^ 



wi 2 ^ wi 

« ^ + f / 

wi 2 1/ 



S(w) 




After projecting the flux J^{w) in the normal direction n = {ux, Uy) (face normal), one can 
compute the Jacobian matrix A„. Its expression in physical variables has the following 
form: 

d{j'{w)-n) I ^ nx riy 

An = — — 7{ = -uun + gHn^ M„ + un^ uriy 

\—VUn + gHny VUrc Un + VUy 
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where Un = uux + vny is the velocity vector projected on n. The Jacobian matrix A„ has 
three distinct eigenvalues: 

\l=Un- C, A2 = Un, A3 = U„ + C, (2.4) 

where c = ^/gH is the speed of gravity waves in the limit of infinite wavelength. This 
quantity plays the same role as the sound speed in compressible fluid mechanics. It is 
now obvious that the system (2.1), (2.2) is strictly hyperbolic provided that H > 0. The 
eigenstructure of the Jacobian matrix A„ is fundamental for constructing the numerical 
flux function (see Section 3.1) and thus, upwinding the discrete solution. 

3. Discretization procedure 

In this study we selected the most natural numerical method for this type of equations. 
Finite volume (FV) methods are a class of discretization schemes that have proven highly 
successful in solving numerically a wide class of systems of conservation laws. These 
systems often come from compressible fluid dynamics. In electromagnetism, for example, 
discontinuous Galerkin methods have proven to be more efficient [CLS04]. When compared 
to other discretization methods such as finite elements or finite differences, the primary 
advantages of FV methods are robustness, applicability on general unstructured meshes, 
and the intrinsic local conservation properties. Hence, with this type of discretization, 
mass, momentum and total energy are conserved exactly, at least in the absence of source 
terms and appropriate boundary conditions. 

In order to solve numerically the system of balance laws (2.1), (2.2) one uses again the 
conservative form of governing equations (2.3). System (2.3) should be provided with an 
initial condition 

w{x,0) = Wo{x), x={x,y)EQ (3.1) 

and appropriate boundary conditions. The implementation of different boundary condi- 
tions will be discussed below (see Section 3.7). 




Figure 2. An example of control volume K with barycenter O. The normal 
pointing from i^' to L is denoted by ukl- 
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3.1. First order scheme. The computational domain i7 C is triangulated into a set 
of non overlapping control volumes that completely cover the domain. Let T denote a 
tesselation of the domain Vt with control volume K such that 

VJKerK = Cl, K ■.= KU dK. 

For two distinct control volumes K and L in T, the intersection is an edge with oriented 
normal ukl or else a vertex. We need to introduce the following notation for the neigh- 
bourhood of K: 

M{K) ■={LeT: aieaiK n L) ^ 0} , 

a set of all control volumes L which share an edge in 2D or a face in 3D with the given 
volume K. In this study, we denote by vol(-) and area(-) the area and length respectively. 

The choice of control volume tesselation is flexible in the FV method. In the present 
study we selected the cell-centered approach (see Figure 3), which means that degrees of 
freedom are associated to cell barycenters. 

/-''T'^----,,,,,,^^^ . Storage location 

. • ■ Control volume 



Figure 3. Illustration for cell-centered finite volume method 

The first steps in FV methods are classical. One starts by integrating equation (2.3) on 
the control volume K shown in Figure 2 and one applies Gauss-Ostrogradsky theorem for 
advective and diffusive fluxes. Then, in each control volume, an integral conservation law 
is imposed: 

^ [ w dVL+ [ T{w) ■ ukl da= [ S{w) dVt (3.2) 
Jk JdK Jk 

Physically an integral conservation law states that the rate of change of the total amount of 

a quantity (for example: mass, momentum, total energy) with density w in a fixed control 

volume K is balanced by the flux J-" of the quantity through the boundary dK and the 

production of this quantity S inside the control volume. 

The next step consists in introducing the control volume cell average for each K 
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After the averaging step, the FV method can be interpreted as producing a system of 
evolution equations for cell averages, since 

' w{x,t) dn = vol{K)- 



dt Jj^ dt 

Godunov was first [God59] to pursue and apply these ideas to the discretization of the gas 
dynamics equations. 

However, the averaging process implies piecewise constant solution representation in each 
control volume with value equal to the cell average. The use of such a representation makes 
the numerical solution multivalued at control volume interfaces. Thereby the calculation 
of the fluxes Jg^{J^{w) ■ ukl) da at these interfaces is ambiguous. A fundamental aspect 
of FV methods is the idea of substituting the true flux at interfaces by a numerical flux 
function 

{T{w) ■ n) ^ ^WK, Wl; Ukl) : x ^ M^ 

a Lipschitz continuous function of the two interface states wk and wl- The key ingredient 
is the choice of the numerical flux function $. In general this function is calculated as an 
exact or even better approximate local solution of the Riemann problem posed at these 
interfaces. In the present study we implemented several numerical fluxes (HLL, HLLC, 
FVCF) described below. 

Any numerical flux is assumed to satisfy the following properties: 

Conservation.: This property ensures that fluxes from adjacent control volumes 
sharing an interface exactly cancel when summed. This is achieved if the numerical 
flux function satisfies the identity 

^{wk^wl^ukl) = -^{wl^wk^ulk)- 

Consistency.: Consistency is obtained when the numerical flux with identical state 
arguments reduces to the true flux of the same state, i.e. 

w; n) = {J-'{w) ■ n){w). 

In the following paragraphs 3.1.1 - 3.1.3 we give several examples of numerical flux 
functions $ which were implemented in the VOLNA code. These choices are justified by 
efficiency, clarity and personal preferences of the authors. However, we do not impose them 
and a final user can easily implement his favourite numerical fiux function. 

3.1.1. FVCF approach. First we describe the scheme called Finite Volumes with Charac- 
teristic Flux (FVCF) and proposed by Ghidaglia et al. in [Ghi95, GKC96, GKCOl]. 
Consider a general system of conservation laws in ID that can be written as follows: 

where ^ G and / : ^ M™. We denote by Aiw) the Jacobian matrix 

and we deal with the case where (3.3) is smoothly hyperbolic, that is to say: for every w 

there exists a smooth basis (ri(w), . . . ,rm{w)) of consisting of eigenvectors of A{w). 
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That is 3Xk{w) G M such that A{w)rk{w) = \k{w)rk{w). It is then possible to construct 
{li{w), . . . , lm{w)) such that ''A{w)lk{w) = Xk{w)lk{w) and lk{w) ■ rp{w) = 5k,p. 

Let M = Ujgz[xj_i/2, a;j+i/2] be a ID mesh. The goal is to discretize (3.3) by a FV 
method. We set Axj = Xj+1/2 — a;j_i/2, At„ = tn+i — tn (we also have IR+ = U„gN[^n, Wi]) 
and 

u;" = ^— / ti;(x,t„)dx, /I+1/2 = ^ / f{w{xj+i/2,t))dt. 
With these notations, we deduce from (3.3) the exaci relation: 

= w^-$^( f-^,,, - fly,) . (3.4) 



J J 



Since the (/j!,_^/2)iGZ cannot be expressed in terms of the (w")jgz, one has to make an 
approximation. In order to keep a compact stencil, it is more efficient to use a three point 
scheme: the physical flux /jYi/2 approximated by a numerical flux 5f"(w", w"^_^). Let us 

show how this flux is constructed here. Since A(w)^ = ^^j^^ we observe that according 
to (3.3) 

J^ + AhJ^^O. (3,5) 

This shows that the flux f{w) is advected hy A{w) like w. The numerical flux g'jiw^ ""^j+i) 
represents the flux at an interface. Using a mean value Aij+1/2 ^^^i^ interface, we 

replace (3.5) by the linearization: 

^^^'^?-)^-. (-) 

It follows that, defining the k-ih. characteristic flux component to be fk{w) = hif^^^i/,) ' 
f{w), one has 

^ + A^(^W^-0. (3.7) 
This linear equation can be solved explicitly: 

fk (W) {X, t) = fk{w) {X - Xk{f^]+l/2) it - tn), tn) • (3.8) 

From this equation it is then natural to introduce the following definition. 

Definition 1. For the conservative system (3.3), at the interface between the two cells 
[xj_i/2,Xj+i/2] and [xj+i/2,Xj+3/2], the characteristic flux g^^ is defined by the following 
formula for A; G {1, . . . , m} : 

we take /i"+i/2 = {^XjiuJ + Axj+iWj^-^) / [Axj + Ax^+i) j 

4(/i^+i/2) ■ <+i) = 4(/i^+i/2) ■ /«) , when Xkif^]^,/2) > , 

^fc(/ii+i/2) • g^'^^i^l = 4(/x^+i/2) ■ > when Xk{fi]+,/2) < , (3.9) 
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when \k{fJ']+i/2) = 0- 

Remark 2. At first glance, the derivation of (3.5) from (3.3) is only valid for continuous 
solutions since A{w) is a non conservative product. In fact, equation (3.5) can he 
justified even in the case of shocks as proved in [Ghi98]. Let us briefly recall here the 
key point. Assuming that the solution undergoes a discontinuity along a family of disjoint 
curves, we can focus on one of these curves that we parameterize by the time variable t. 
Hence, locally, on each side of this curve, w{x,t) is smooth and jumps across the curve 
X = ^{t). The Rankine-Hugoniot condition implies that f{w{x,t)) — a{t)w{x,t) , where 
a if) = is smooth across the discontinuity curve and therefore A{w)^^^^ can be 

defined as A{w)^-^ = A{w) ^^^^"gj;""'^ + ■ 

Proposition 1. Formula (3.9) can be written as follows: g^^'"'{w'j, = g^^{^'^] w", w^^-^ 

where 



\ . i \ , ( ,,\—f\ V 



f{v) + /H , 



\k(p)<0 Afc(M)=0 

+ E ilM-fiv))r,ifi). (3.10) 

Afc(M)>0 

Proof. This comes from the useful identity vahd for all vectors $ and fi in M™: 

k=m 

$ = ^^(/A;(yu) ■ ^)i"kilJ')- We also observe that (3.10) can be written under the following 

k=l 

condensed form: 

g {fi;v,w) = U{fi;v,w) , (3.11) 

where f/(/i; v, w) is the sign of the matrix A{fi) which is defined by 

m 

sign(A(/i))$ = J]sign(Afc)(/fc(^) • $)rfc(/i). 

k=l 

The form (3.11) refers to a numerical flux leading to a flux scheme [Ghi98]. ■ 

Remark 3. Let us discuss the relation, in the conservative case, between the characteristic 
numerical flux g^^ and the numerical flux leading to Roe's scheme [RoeSl]. The latter 
scheme relies on an algebraic property of the continuous flux f{w) which is as follows. It 
is assumed that for all admissible states v and w, there exists a m x m matrix A^^^{v, w) 
such that f{v) — f{w) = A^^^{v,w){v — w) (Roe's identity). Then the numerical flux 
leading to Roe's scheme is given by: 
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But using Roe's identity, we obtain that 



/H - f{v) 



2 



(3.13) 



which is of the form (3.11): Roe's scheme is also a flux scheme. The characteristic flux 
proposed in this paper is more versatile than Roe's scheme in the sense that it does not rely 
on an algebraic property of the flux. Hence for complex systems (like those encountered in 
the context of two phase flows) this scheme is an efficient generalization of Roe's scheme. 
Moreover, as we shall see below, this scheme has a natural generalization to arbitrary non 
conservative systems. Finally, the fact that the numerical flux is a linear combination of 
the two fluxes induces a quite weak dependence on the state fi which appears in formula 
(3.10), see [CGOO]. 

3.1.2. HLL numerical flux. Now we present another approximate Riemann solver which 
was proposed by Harten, Lax and van Leer [HLvL83]. Nowadays this method is known as 
the HLL scheme. While the exact solution to the Riemann problem contains a large amount 
of detail, the HLL solver assumes fewer intermediate waves. The simplified Riemann fan 
is illustrated on Figure 4. It consists of two waves separating three constant states. 



sl t 








X 



Figure 4. Approximate Riemann fan corresponding to the HLL scheme. 



Consider the following Riemann problem: 




X < 0, 
X > 0. 



(3.14) 



Wr 




Figure 5. Two states wl and wr connected by Rankine-Hugoniot curves 
represented in the phase space. 
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The intermediate state in the approximate Riemann fan will be denoted by w* and 
two shock wave speeds are denoted by sl and respectively (see Figures 4 and 5 for 
illustration). In order to determine the unknown intermediate state, we write the Rankine- 
Hugoniot conditions twice: 

sl{w*-wl) =F*-Fl, 
SR{wn-w*) =Fr-F\ 

where -Fl,_r := F{wl,r)- It is straightforward to find the solution to this system: 

SrWr - SlWl - {Fr - Fl) 



w 



sr - Sl 




P , / * x SrFl - SlFr + SlSr{wr - Wl) ,^ ^ 

F =Fl + sl{w -wl) = . (3.15) 

Sr - Sl 

Now we have all the elements to define the numerical flux of the HLL scheme: 

^hll{wl,Wr) : = 

During the presentation of the HLL scheme we missed one important point: how to 
estimate the wave speeds sl and srI The answer is crucial for the overall performance of 
the scheme. With appropriate choices for the wave speeds sl and sr, the HLL scheme pos- 
sesses nice numerical properties. Namely, it satisfies an entropy inequality [Dav88], resolves 
isolated shocks exactly [HLvL83] and preserves positivity [EMRS91]. In our code we im- 
plemented the following choice for sl and sr which is motivated by analytical expressions 
for the Jacobian eigenvalues (2.4): 

Sl = mm{uL - Cl,u* - c*), Sr = min{u* + c*,Ur + Cr), 



where cl,r '■= ^/gHLji is the gravity wave speed for the left and right states and 
1 11 

U* = -{ul + Ur) +Cl- Cr, C* = -{cl + Cr) - -{ur - Ul). 

Numerical experiments show that this approximate Riemann solver is very robust with the 
above choice for the wave speeds [CIM+OO, ZCIM02]. One can show that the HLL scheme 
belongs to the class of flux schemes. Recall that a FV scheme is called a flux scheme if its 
numerical flux can be written in the following form: 



2 ^ 2 

where U{wl,wr) is some matrix. The robustness of the HLL scheme can be explained by 
this nice property. 

However, the HLL scheme has one important shortcoming: it cannot resolve isolated 
contact discontinuities. In the next section 3.1.3 we present another scheme which was 
designed to remedy this problem. 
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3.1.3. HLLC flux. The HLL scheme presented briefly in the previous section was later 
improved by Toro, Spruce and Speares [TSS94]. Their modification concerns essentially 
the structure of the Riemann fan which is depicted on Figure 6. Namely, they introduced a 
contact discontinuity between two shock waves of the HLL scheme. That is why the novel 
scheme was called the HLLC scheme [FT95]. 



sl 


t 


s* 












//-^ Wr 







Figure 6. Approximate Riemann fan corresponding to the HLLC scheme. 



Here we do not provide details on the derivation of the HLLC scheme and refer to the 
original articles and others which can fill this gap [BCCC97, KCY07]. 

We consider the same Riemann problem (3.14). In the HLLC approximation, the solution 
to this Riemann problem consists of three waves with speeds sl, s* and sr separating four 
constant states wl, w^, Wr and wr. Wave speeds sl,r are estimated as in previous section 
3.1.2, while s* is given by the formula 

SlHr{ur - Sr) - SrHl{ul - Sl) 

S = . 

Hr{ur - Sr) - Hl{ul - Sl) 

The intermediate states w^ = j^, {Hu)\ j^, {Hv)\ j^) are computed as follows: 

TT* Sl,B~Ul.R TT 

^L,R - SLR-s* ^L,R, 

{hu)Ir = '-'^{hu)l,r + fg!,H ^"""i:;:^^^if'"""\ 

sl,r-ul,r ( u„,\ I g Tj2 {'isL,R-s'' -ul,r){s* -ul^r) 



iHv)lR=''^^iHv)L,R+IH[ 



2^^L,R (si,_R-s*)3 



Finally, the numerical flux of the HLLC scheme is defined as 
^hllc{wl,Wr) : = 



Fl, Sl > 0, 

F*:=FL + SLiwl-WL), Sl<0<s*, 

F^:= Fr + SR{w*ji-WR), s*<0<Sr, 

Fr, Sr < 0. 



3.2. Semidiscrete scheme. Introducing the cell averages wk and numerical fluxes into 
(3.2) yields for the integral conservation law 

dwK area(L n K) ^ , ^ , 1 f . ,^ 

^+ E \p^HwK,WL-,nKL) = —r7JT, S{w)dn. 

at ^ vo1(a) vo1(a) 



LeN{K) 



We denote by Sk the approximation of the quantity ^^^^^ S{w) dfl. The source term 
discretization is discussed in Section 3.4. Thus, the following system of ordinary differential 
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equations (ODE) is called a semi-discrete FV method: 

dwK area(L n K) ^ , ^ . „ , ^ , , 

Yl voim ^ i^K,WL;nKL)=SK, yKeT. (3.16) 

The initial condition for this system is given by projecting (3.1) onto the space of piecewise 
constant functions 

1 f 



vol K J 

This system of ODE should also be discretized. There is a variety of explicit and implicit 
time integration methods. Let denote a numerical approximation of the cell average 
solution in the control volume K at time t" = nAt. The simplest time integration method 
is the forward Euler scheme 

dwK ^ wl+^ - wl 
dt At 

When applied to (3.16) it produces the fully-discrete FV scheme: 

— \ area(L n K) ^ , „ „ ^ , ^„ , ^ 

+ E voim H<,wl;nKL) = S^^, VKeT. (3.17) 

The time discretization used in this study is detailed in Section 3.5. 

3.3. Run-up algorithm. As already pointed out above, the NSWE are strictly hyperbolic 
if if > 0, i.e. when some water is present. The shoreline position is given by the implicit 
relation H{x,t) = 0. At these locations the system loses its strict hyperbolicity. Finally, 
in dry regions, H < 0, the system is non-hyperbolic, i.e. ill-posed. All these facts mean 
that there are some major theoretical difficulties in considering the inundation problem. 

Very often some ad-hoc artificial techniques are implemented to circumvent run-up and 
run-down problems ("slot technique" of Madsen et al. [MSS97], algorithm of Hibberd- 
Peregrine [HP79], use of coordinate transformations [OHK97] and so on) - see also [MBFS07]. 

The shoreline boundary conditions have a very simple analytical form: 

dx 

H{xs{t),t) = 0, -^=u{xs{t),t), 

where Xs{t) is the shoreline position. 

The algorithm proposed by Brocchini et al. [BBMAOl] was chosen for the VOLNA code. 
It is based on the shoreline Riemann problem shown in Figure 7 [Sto57]: 

{dw I dFjw) ^ p ( dw , dF{w) ^ ^ 

dt 9x ' 9t dx ' 

The main idea to solve the shoreline Riemann problem is to pass to the limit wl — or 
wr — )■ in the solution to the classical Riemann problem (3.14). Technical details can 
be found in [BBMAOl]. However, we do not need to know the complete solution. It is 
sufficient to extract the wave propagation speeds at the shoreline (see Figure 8). These 
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Figure 7. Shoreline left and right Riemann problem. 



analytically determined speeds are imposed in an approximate Riemann solver when a 
wet/dry transition is detected. 

Consider two control volumes K and L which share a common face. We must find the 
numerical fiux ^{wK.WL.nKLj across this face. Let us summarize the key points of the 
method: 

Wet/wet interface:: If Hl > and > 0, we apply in the usual way an approx- 
imate Riemann solver which gives the numerical fiux 

Dry/dry interface:: If Hi = Hr = 0, we just return the zero fiux $ = since there 
is no fiow between two dry cells. 

Wet/dry interface:: If Hji = and > 0, we have a situation corresponding to 
the left shoreline Riemann problem. Its solution yields the following choice of the 
wave speeds: 

sl ■= {ul ■ ukl) - cl, sr := {ul ■ ukl) + 2cl. 

Then we apply the HLL or the HLLC scheme with the above values of sl and sr 
(see Figure 8 for illustration). 
Dry/wet interface:: If Hr > and Hi = 0, we have a situation symmetric to the 
previous case. One must solve the right shoreline Riemann problem. It provides 
the following speeds: 

Sl ■= {ur ■ Ukl) - "^cr, sr := {ur ■ ukl) + cr. 
Here again, the HLL or the HLLC scheme is applied. 

We would like to underline the simplicity of this approach. In fact, there is no special 
treatment for the interface. This algorithm is run uniformly in the whole computational 
domain leading to an easy and robust implementation. We validate this method in sections 
4.2 - 4.4. 

3.4. Source terms discretization. In this section we discuss some issues related to the 
source term discretization and we explain a technique to remedy them. 
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Figure 8. Shoreline Riemann problem (left) and wave propagation speeds. 

Source terms of the form gHVh arise in the horizontal momentum conservation equation 
(2.2). Obviously, this term is equal to zero when the bottom is flat {h = const). However, it 
is not the case in real world applications. The magnitude of this term is proportional to the 
bed slope and may take large values when abrupt changes are present in the bathymetry. 

Another profound property of NSWE is that the system (2.1), (2.2) admits non-trivial 
steady states. They can be determined from the following steady equations: 

r V ■ (Hu) = 0, 

\ V ■ {Hu®u + fif2) = gHVh 

It is not so trivial to find analytical solutions to these equations. However, an ideal numer- 
ical scheme should preserve them. Recall that for ID flows it is possible to describe the 
whole family of steady states and this information can be used to design efficient source 
term discretizations [LR98, VC99]. Since most applications require a 2D solver, we address 
this problem directly in 2D. 

As explained above, it seems to be extremely difficult to construct a scheme which 
preserves exactly all steady state solutions. Thus, we have to simplify the problem. We will 
focus our attention on a simple class of steady solutions which are called in the literature 
"lake at rest": 

M = 0, 1] := H — h = const. (3.18) 
The last relations can be expressed in discrete variables: 

uk = ul = 0, Hk — hx = Hl ~ h-L = const. (3.19) 

We briefly present the method chosen for our code and developed in [ABB"'"04, AB05]. It 
is based on the idea of the interface hydrostatic reconstruction. 

The well-balanced algorithm takes as input the vector of conservative variables {wA'jA'eT; 
bathymetry data {hK}KeT is composed of the following steps: 

• Assume that the control volumes K and L share a common face K (1 L. In this 
case, the interface bathymetry is deflned as := m.m{hK, hi). This step is done 
only once at the initialization stage. 
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• The hydrostatic reconstructed interface water depth is given by 

H^j^ = {Hk — hx + h*j^L) + ^ where = ma.x{z, 0). 

From the dicrete interpretation (3.19) of the well-balanced condition (3.18), we 
define a new vector of the interface conservative variables: 

- (4^5.) ■ <^-'»* 

• From the balance of hydrostatic forces V(|-ff^) = gHVh, the adapted discretiza- 
tion of the source terms is introduced: 

• The well-balanced scheme is obtained by replacing cell-centered values wk by new 
interface values (3.20): 

w'^^-w]^ ^ area(Lnir) ^ 

At + vol(K) ^(^KM;nKL) = 

It can be proven [AB05] that the hydrostatic reconstruction strategy preserves the "lake 
at rest" solutions and ensures the positivity property. We describe here only the first-order 
algorithm for the sake of simplicity. The extension to second order can be found in the 
original papers and in [Aud04]. 

3.5. Time discretization. In the previous sections we considered the spatial discretiza- 
tion procedure with a FV scheme. It is a common practice in solving time- dependent PDEs 
to first discretize the spatial variables. This approach is called method of lines: 

wt + dj{w) = S{w) ^wt = C{w) (3.21) 

In order to obtain a fully discrete scheme, we must discretize the time evolution operator. 
In the present work we chose the so-called Strong Stability-Preserving (SSP) time dis- 
cretization methods described in [Shu88, GSTOl, SR02]. Historically these methods were 
called Total Variation Diminishing (TVD) time discretizations. 

The main idea behind SSP methods is to assume that the first order forward Euler 
method is strongly stable (see the definition below) under a certain norm for the method 
of lines ODE (3.21). Then, we try to find a higher order scheme. Usually the relevant 
norm is the total variation^ norm: 

TV(^") :=5^|<-<-i| 

3 

and TVD discretizations have the property TV(w"'''"^) < TV{w^). 



The notion of total variation is used essentially for ID discrete solutions. 



^(0) 


= W", 












fc=0 
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Remark 4. Special approaches are needed for hyperbolic PDEs since they contain discon- 
tinuous solutions and the usual linear stability analysis is inadequate. Thus a stronger 
measure of stability is usually required: 

Definition 2. A sequence {w"} is said to be strongly stable in a given norm ||-|| provided 
that I I < I I for all n > 0. 

A general m-stage Runge-Kutta method for (3.21) can be written in the form 

(3.22) 

'=) + At/3i,fci:(ti;«)), ai,k>0, i = l,...,m, (3.23) 

(3.24) 

In [S088] the following result is proved 

Theorem 1. // the forward Euler method is strongly stable under the CFL restriction 
At < AtpE 

+ At£(w")|| < 

then the Runge-Kutta method (3.22) - (3.24) with Pi^k > zs SSP, Ww'^+^W < Ww'^W, 
provided the following CFL restriction is fulfilled: 

At<cAtFE, c = min— 

Pi,k 

Here we give a few examples of SSP schemes which are commonly used in applications 
(optimality is in the sense of CFL condition): 

• Optimal second order two-stage SSP-RK(2,2) scheme with CFL = 1: 



Zi Zi 

Optimal third order three-stage SSP-RK(3,3) scheme with CFL = 1: 

^(1) = + Ati:(w(")), 

4 4 4 ^ ^' 

o o o 
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• Third order four-stage SSP-RK(3,4) scheme with CFL = 2: 
^(1) = + iAt£(w(")), 

3 3 6 

^(n+l) ^ + iAt£(w(3))_ 

The hnear absolute stabihty region for the RK and SSP-RK schemes is the same. However 
the nonhnear absolute stability regions are quite different [CP92]. 

We tested these different schemes in our numerical code and decided to adopt SSP- 
RK(3,4) due to its accuracy and wide stability region. In our opinion this scheme represents 
a very good trade-off between precision and robustness. 

3.6. Second order extension. If we analyze the above scheme, we understand that in 
fact, we have only one degree of freedom per data storage location. Hence, it seems that 
we can expect to be first order accurate at most. In the numerical community first order 
schemes are generally considered to be too inaccurate for most quantitative calculations. Of 
course, we can always make the mesh spacing extremely small but it cannot be a solution 
since it makes the scheme inefficient. From the theoretical point of view the situation 
is even worse since an 0{h2) Li-norm error bound for the monotone and E-flux schemes 
[Osh84] is known to be sharp [Pet91], although an 0{h) solution error is routinely observed 
in numerical experiments. On the other hand, Godunov has shown [God59] that all linear 
schemes that preserve solution monotonicity are at most first order accurate. This rather 
negative result suggests that a higher order accurate scheme has to be essentially nonlinear 
in order to attain simultaneously a monotone resolution of discontinuities and high order 
accuracy in continuous regions. 

A significant breakthrough in the generalization of FV methods to higher order accu- 
racy is due to N.E. Kolgan [Kol72, Kol75] and van Leer [vL79]. They proposed a kind of 
post-treatment procedure currently known as solution reconstruction or MUSCL (Mono- 
tone Upstream-centered Scheme for Conservation Laws) scheme. In the above papers the 
authors used linear reconstruction (it will be chosen in this study as well) but this method 
has already been extended to quadratic approximations in each cell [BF90]. 

3.6.1. Historical remark. In general, authors of numerical articles which use the MUSCL 
scheme often cite the paper by van Leer [vL79]. It is commonly believed in the scientific 
community that B. van Leer was first to propose the gradient reconstruction and slope 
limiting ideas. Because of unfortunate political reasons, N.E. Kolgan's work [Kol72, Kol75] 
remained unknown for a long time. We would like to underline the fact that the first 
publication of Kolgan came out seven years before van Leer's paper. Van Leer seems to be 
aware of this situation since in his recent review paper [vL06] one can find "A historical 
injustice" section: 
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"It has been pointed out to me by Dr. Vladimir Sabelnikov, formerly of 
TsAGI, the Central Aerodynamical National Laboratory near Moscow, that 
a scheme closely resembling MUSCL (including limiting) was developed in 
this laboratory by V. P. Kolgan (1972). Kolgan died young; his work ap- 
parently received little notice outside TsAGI." 

3.6.2. TVD and MUSCL schemes. There is a property of scalar nonlinear conservation 
laws, which was probably observed for the first time by P. Lax [Lax 73]: The total increasing 
and decreasing variations of a differentiable solution between any pair of characteristics 
are conserved. In the presence of shock waves, information is lost and the total variation 
decreases. For compactly supported or periodic solutions, one can establish the following 
inequality 

+ 00 +00 

j \dw{x,t2)\< j \dw{x,ti)\, t2>ti. (3.25) 

— oo — oo 

This motivated Harten [IIar83] to introduce the notion of discrete total variation of nu- 
merical solution Wh '■= {wj} 

TV{wh) ■■= ^ \wj+i -Wj\, 
j 

and the discrete counterpart to (3.25) 

TV«+i) < Tr«). 

If this property is fulfilled, then a FV scheme is said to be total variation diminishing 
(TVD). The following theorem was proved in [IIar83]: 

Theorem 2. (i) Monotone schemes are TVD; (ii) TVD schemes are monotonicity pre- 
serving, i.e. the number of solution extrema is preserved in time. 

Remark 5. From the mathematical point of view it would be more correct to say "the total 
variation non-increasing (TVNI) scheme" but the "wrong" term TVD is generally accepted 
in the scientific literature. 

In one space dimension the construction of TVD schemes is not a problem anymore. Let 
us recall that in this study we are rather interested in two space dimensions (or even three 
in future work). In these cases the situation is considerably more complicated. Even if 
we consider the simplest case of structured cartesian meshes and apply a ID TVD scheme 
on a dimension-by-dimension basis, a result of Goodman and Leveque shows [GV85] that 
TVD schemes in two or more space dimensions are only first order accurate. Motivated by 
this negative result, weaker conditions yielding solution monotonicity preservation should 
be developed. 

In this article we describe the construction and practical implementation of a second- 
order nonlinear scheme on unstructured (possibly highly distorted) meshes. The main 
idea is to find our solution as a piecewise affine function on each cell. This kind of linear 
reconstruction operators on simplicial control volumes often exploit the fact that the cell 
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average is also a pointwise value of any valid (conservative) linear reconstruction evaluated 
at the center of a simplex. This reduces the reconstruction problem to that of gradient 
estimation given cell averaged data. In this case, we express the reconstruction in the form 

wxix) = wk + C^w)k ■ (x - xo), KeT, (3.26) 

where wk is the cell averaged value given by the FV method, {'Vw)k is the solution gradient 
estimate (to be determined) on the cell K, x & K and the point Xq is chosen to be the 
center for the simplex K. 

It is important to note that with this type of representation (3.26) we remain absolutely 
conservative, i.e. 

— ^— I wk{x) dVt = Wk 
■vo\[K) Jk 

due to the choice of the point xq. This point is crucial for FVs because of intrinsic conser- 
vative properties of this method. 

In the next sections we describe briefly two common techniques: Green-Gauss integration 
and least squares methods for solution gradient estimation on each cell. There are other 
available techniques. We can mention here an implicit gradient reconstruction method 
proposed in [MG96] and reused later in [AMS04]. We decided not to implement this 
approach in our research code since this procedure is computationally expensive'^. 




Figure 9. Illustration for Green-Gauss gradient reconstruction. Control 
volume K with barycenter O and exterior normal n. 

3.6.3. Green-Gauss gradient reconstruction. This gradient reconstruction technique can be 
easily implemented on simplicial meshes. It is based on two very simple ideas: the mean 
value approximation and Green-Gauss-Ostrogradsky formula. 

Consider a control volume K with barycenter O. The exterior normal to an edge e G dK 
is denoted by fie. This configuration is depicted on Figure 9. In order to estimate the 



In order to reconstruct the solution gradient we have to solve a linear system of equations. Recall that 
the gradient is estimated at each time step on each control volume. This factor slows down considerably 
explicit time discretizations. 
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solution gradient on K (or in other words, to estimate its value at the center O) we make 
the following mean value approximation 

(V..k=(V.OIo^;;;i^Xv„,da 

and apply Green-Gauss-Ostrogradsky formula 

{Vw)k = — fTTTT / w®nda = } [ w ^ He da = 

vol{K) yol{K) Je 



eGdK ^ ' 

where w|g/2 denote the solution value at the face (or edge in 2D) centroid. The face value 
needed to compute the reconstruction gradient can be obtained from a weighted average 
of the values at the vertices on the face [HC89]. In 2D it simply becomes 

= ^ • 

This approximation yields the following formula for gradient estimation: 

^-^ vol A 2 

The gradient calculation is exact whenever the numerical solution varies linearly over the 
support of the reconstruction. 

This procedure requires the knowledge of the solution values at the mesh nodes {iVj}. 
Since a cell centered FV scheme provides data located at cell centers, an interpolation 
technique is needed. The quality of Green-Gauss gradient reconstruction greatly depends 
on the chosen interpolation method. The method chosen here is explained in Section 3.6.6. 




Figure 10. Illustration for least-squares gradient reconstruction. A trian- 
gle control volume with three adjacent neighbors is depicted. 
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3.6.4. Least-squares gradient reconstruction method. In this section we consider a triangle'^ 
control volume K with three adjacent neighbors Ti, T2 and T3. Their barycenters are 
denoted by 0(xo), Oi(xi), 02(x2) and 03(x3) respectively. In the following we denote by 
Wi the solution value at the centers Of. 

Wi:=w{xi), Wo:=w{xo). 

Our purpose here is to estimate Vw = [d^w, dyw) on the cell K. Using Taylor formula, 
we can write down the three following relations: 

w,-wo = {yw)K-{x,-xo) + 0{h^), z = l,2,3. (3.27) 

If we drop higher order terms 0{h'^), these relations can be viewed as a linear system of 
three equations for two unknowns^ {dxW,dyw). This situation is due to the fact that the 
number of edges incident to a simplex mesh in is greater or equal (in this case see 
Remark 6) to d thereby producing linear constraint equations (3.27) which will be solved 
analytically here in a least squares sense. 

First of all, each constraint (3.27) is multiplied by a weight Ui G (0, 1) which will be 
chosen below to account for distorted meshes. In matrix form our non-square system 
becomes 

/uiAxi wiAyA /wi(wi 

UJ2AX2 W2Ay2 C^w)k = ^2(^2 

ywaAxs Us Ays J ^^3(^3 

where Axj = Xi — xq, Ayi = yi — yo- For further developments it is convenient to rewrite 
our constraints in abstract form 

[Li, L2] ■ {Vw)k = f. (3.28) 

We use a normal equation technique in order to solve symbolically this abstract form in a 
least squares sense. Multiplying on the left both sides of (3.28) by [L1L2]* yields 

GiVu.U.i, G ^ ^((|;^;j l^-g) (3.29) 

where G is the Gram matrix of vectors I Lt,L2\ and b = \^^^ ■ The so-called normal 
equation (3.29) is easily solved by Cramer's rule to give the following result 

(yw)K - — w- 




klk2 — ^12 V^ll(-^2 ■ /) — ^12(-^l ■ /). 

The form of this solution suggests that the least squares linear reconstruction can be 
efficiently computed without the need for storing a non-square matrix. 

Now we discuss the choice of weight coefficients The basic idea is to attribute 

bigger weights to cells barycenters closer to the node iV under consideration. One of the 



^Generalization to other simplicial control volumes is straightforward. 

^This simple estimation is done for the scalar case only w = (w). For more general vector problems the 
numbers of equations and unknowns must be changed depending on the dimension of vector w. 
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possible choices consists in taking a harmonic mean of respective distances rj = — xatI 
This purely metric argument takes the following mathematical form: 

-k 



E"^ 1 1 ~* ~* 1 1 — fc ' 
7=1 iFj ~ ^n\\ 



where k in practice is taken to be one or two (in our code we choose k = 1). 

Remark 6. When a triangle shares an edge with the boundary dVt (see Figure 12 for 
illustration ), the gradient reconstruction procedure becomes even simpler, since the number 
of constraints is equal to d and the linear system (3.27) becomes completely determined: 

Wi-wo = (Vw)if ■ (fi - xo) + i = l,2. 

In component form it reads 

xi -Xq yi- yo\ , . _ (wi-Wq 



X2-XQ 1/2-1/0/ - Wo 

The unique solution to this linear system is given again by Cramer's rule 

(l/2 - yo){.wi - wq) - {yi - yo)iw2 - wq) 

(Xi - Xo){w2 - Wo) - ix2 - Xo){Wi - Wq) 



K 



{xi - xo){y2 - yo) - {x2 - xo){yi - yo) 

3.6.5. Slope limiter. The idea of incorporating limiter functions to obtain non-oscillatory 
resolution of discontinuities and steep gradients goes back to Boris and Book [BB73]. When 
the limiter is identically equal to 1, we have the unlimited form of the linear interpolation. 
In the ID case one can easily find in the literature about 15 different limiter functions such 
as CHARM, minmod, superbee, van Albada and many others. On unstructured meshes the 
situation is quite different. In the present study we decided to choose the Barth-Jespersen 
limiter proposed in [BJ89]. Here we do not discuss its construction and properties but just 
give the final formula. We need to introduce the following notation 

wT"" := min wl, w^""" := max wl ■ 

L€M{K) L€M{K) 

The limited version of (3.26) is given by the following modified reconstruction operator 

Wxix) = Wk + aKiyw)K ■ {x - Xq), K , 

where it is assumed that ax G [0, 1]. Obviously, the choice ax = corresponds to the first 
order scheme while ax = 1 is the unlimited form. Barth and Jespersen [BJ89] propose the 
following choice of ax- 

1 otherwise, 

where Xf denotes the face / centroid. 

Although this limiter function does not fulfill all the requirements of FV maximum prin- 
ciple on unstructured meshes [BO04], it can be shown that it yields FV schemes possessing 
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a global extremum diminishing property. Also this limiter produces the least amount of 
slope reduction which can be advantageous for accuracy. Note that in practice minor 
modifications are required to prevent near zero division for almost constant solution data. 

3.6.6. Solution interpolation to mesh nodes. We have seen above that several gradient re- 
construction procedures (in particular gradient estimation on the faces) require the knowl- 
edge of the solution at mesh nodes (or vertices). This information is not directly given by 
the FV method since we chose the cell-centered approach. 




Figure 11. Triangles with their barycenters Oj sharing the same vertex A^. 



Let us consider a node N{xn, Un) of the tesselation T and a control volume Ki with 
barycenter Oi{xi,yi) having this node as a vertex (see Figure 11 for illustration). The 
MUSCL procedure provides a solution gradient on each cell. Thus, using the Taylor formula 
or, equivalently, the representation (3.26) we can estimate the solution value at the node 
N 

wn = wk, + (yw)K, ■ {xn - Xi). (3.30) 

The problem is that we will have d{N) different values of the solution in the same point 
depending on the control volume under consideration. Here d{N) is the degree of vertex 
in the sense of graph theory. One of the possible ways to overcome this contradiction is 
averaging. One interesting technique was proposed in [HC89], further improved in [KMC03] 
and slightly modified by us. The algorithm implemented in our code is briefly described 
here. 

First of all, let us look for the vertex value iDat as a weighted sum of the values 
computed by formula (3.30) from each surrounding cell 

The weighting factors {uJiYi^ made to satisfy the condition of zero pseudo-Laplacian 

d{N) d{N) 

L{Xn) = (^iiXi - Xn), L{yn) = ^ UJiiVi - Vn) ■ (3.31) 

i=l i=l 
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These conditions have a very simple interpretation. They are imposed so that the method 
be exact for affine data over the stenciL 

As in the original formulation by Holmes and Connell [HC89], the weighting factor uji is 
written as 

uji = 1 + Aui . 

The weights {ui} are determined by solving an optimization problem in which the cost- 
function to be minimized is defined as 

-^(r^Au;,)' ^min (3.32) 

i=l 

with two constraints given by (3.31). It should be noted that the cost function is slightly 
different from the original formulation. The difference lies in the factor of 

which was introduced in [KMC03]. This modification effectively allows larger values of 
weight Aui for those cells closer to the node in question. 

Employing the method of Lagrange multipliers, the original optimization problem, which 
was to minimize the cost function given by (3.32) with the constraints (3.31), is equivalent 
to minimizing the function C defined by 



£ = ^ {riAuiY - A ^ Ui{xi - x^) - ^ uji{yi - ?/„) min 

i=l i=l 
X{Xi - Xn) + fiiVi - Vn) 



i=l i=l i=l 

which leads to 

Aui 



j,2 



The two Lagrangian multipliers, A and /i, are obtained from 

■1 ^ylxy ^xlyy fx^xy ^y^xx 

II - J2 ' ^^11 - P ' 

^xx-^yy J-xy ^xx-^yy ^ xy 



where 



d{N) d{N) 

^{Xi-Xn), ry=^{yi-yn) 



i=l i=l 
[Xj- XnY _Sr^ \Vi- Vn) r _ [Xj - XnjjVi - Vn) 

P ' yy p ' ""y p 

i=i « i=i « i=i « 

The last step consists in renormalizing the weights {wjj^f^^^ to the range [0, 1]. 

Remark 7. The above algorithm is not computationally expensive since the weights {uJi}'!^'' 
only depend on the tesselation T geometry. It means that they can be computed and stored 
before the main loop in time and reused during later computations. 
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Figure 12. Control volume sharing a face with boundary dQ. 



3.7. Implementation of boundary conditions. So far we have not discussed the im- 
plementation of boundary conditions. The flavor of the treatment of boundary conditions 
for hyperbolic systems is given here and we refer to [GP05] for a general discussion. This 
is a very important topic since they actually determine the solution. Let us consider the 
space discretization of the system (3.3) by a cell centered FV method. For instance for 
the time explicit discretization we have the scheme (3.17). Of course this formula is not 
valid when K meets the boundary of Q (see Figurel2 for illustration). When this occurs, 
we must find the numerical flux ^{v^, K,dfl). In practice, this flux is not given by the 
physical boundary conditions and moreover, in general, (3.3) is an ill-posed problem if we 
try to impose either w or F{w) ■ n on dVt. This can be understood in a simple way by using 
the following linearization of this system: 

where n represents the direction of the external normal on K r\ dQ, is the advection 
matrix: 

_ {dF{w) -n) 

Aa — 'Qw I«"=m' (3.34) 

and w_ is the state around which the linearization is performed. When (3.3) is hyperbolic, 
the matrix is diagonalizable on M and by a change of coordinates, this system becomes 
an uncoupled set of m advection equations: 

^ + Afc^ = 0, fc = l,...,m. (3.35) 

Here the are the eigenvalues of A^ and according to their sign, waves are going either 
into the domain (A^ < 0) or out of the domain (A^ > 0). Hence we expect that it is only 
possible to impose p conditions on ii" fl (9fi where p = tl{A;G{l,...,m} such that A^ < 0}. 
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Let US consider now a control volume K which meets the boundary dVL. We take w_ = 
and write the previous linearization. We denote by x the coordinate along the outer normal 
so that (3.33) reads: 

dw . dw 

which happens to be the linearization of the ID [i.e. when nd = 1) system. First we label 
the eigenvalues Xk{w) of by increasing order: 

Xi{w) < < . . . < Xpiw) < < Xp+iiw) ...< A„(w) . (3.37) 

(i) The case p = 0. In this case information comes from inside Q and therefore we 
take: 

$(w^, K, an) = F{wl) ■ UK . (3.38) 

In the Computational Fluid Dynamics (CFD) literature this is known as the "su- 
personic outflow" case. 

(ii) The case p = m. In this case information comes from outside Q and therefore we 
take: 

^wl,K,dn) = ^gi,en, (3.39) 

where ^given are the given physical boundary conditions. In the CFD literature 
this is known as the "supersonic inflow" case. 

(iii) The case 1 < p < m — 1. As already discussed, we need p scalar information coming 
from outside of Q. Hence we assume that we have on physical ground p relations 
on the boundary: 

giiw) = 0, 1 = 1,..., p. (3.40) 

Remark 8. The notation gi{w) = means that we have a relation between the components 
ofw. However, in general, the function gi is not given explicitly in terms ofw. For example 
gi{w) could be the pressure which is not, in general, one of the components ofw. 

Since we have to determine the m components of ^{w]^, K,dQ), we need m — p 
supplementary scalar conditions. Let us write them as 

hi{w) = 0, I =p + l,...,m. (3.41) 

In general (3.40) are referred to as "physical boundary conditions" while (3.41) are 
referred to as "numerical boundary conditions" . 
Then we take: 

^{w'^^K^dn) = F{w)-nK, (3.42) 
where w is solution to (3.40)-(3.41) (see however Remark 11 and (3.48)). 

Remark 9. The system (3.40)-(3.41) for the m unknowns w&Gisamxm nonlinear 
system of equations. Its solvability is given by Theorem 3. 
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Let US first discuss the numerical boundary conditions (3.41). By analogy with what we 
did on an interface between two control volumes K and L, we take (recall that w_ = w^): 

h{w) ■ {F{w) ■ uk) = lk{w) ■ {F{wl) -hk), k=p + l,...,m. (3.43) 

In other words, we set hkiw) = hiw'^) ■ {F{w) ■ uk) — hiw^) ■ {F{w^) ■ uk)- We have 
denoted by {h{w), . . . , lm{w)) a set of left eigenvectors of A^: ^AJkiw) = ^kh{w) and by 
{T^iiw.): ■ ■ ■ j^mim)) a set of right eigenvectors of A^: A^rk{w) = Xk^kiw)- Moreover the 
following normalization is taken: lk{w) ■ fp{w) = Sk,p- 

According to [GP05] we have the following result on the solvability of (3.40)-(3.41). 

Theorem 3. In the case 1 < p < m — 1, assume that Xp+i{w) > 0, and 

With the choice (3.43) the nonlinear system (3.40)-(3.41) has one and only one solution 
V, for V — w and gi{w) sufficiently small. 

Remark 10. In this result we exclude the case where the boundary is characteristic i.e. 
the case where one of the Xk is equal to 0. This case cannot be dealt with at this level 
of generality. On the other hand, wall boundary conditions belong to this category. They 
can be discussed and handled directly on the physical system under consideration. In this 
section we show how to do it for the NSWE equations (see Paragraph 3.7.1). Moreover, the 
treatment of wall boundaries of compressible Euler equations and some two-phase systems 
[DDGOSb, DDGIO, DDGOSa, Dut07] can be done in a similar way. 

Remark 11. In practice, (3.40)-(3.41) are written in a parametric way. We have a set of 
m physical variables w (e.g. pressure, densities, velocities,. . .) and we look for w satisfying: 

gi{w)=^, / = l,...,p, (3.45) 

lk{w) ■ $ = lk{w) ■ iF{wl) ■ Uk) , (3.46) 
$ = F{w) ■ Uk , (3.47) 

and then we take: 

^wl,K,dQ) = <^. (3.48) 
The system (3.45)-(3.46)-(3.47) is then solved by Newton's method. 

3.7.1. Impermeable boundary. Consider the case of a rigid wall boundary 

u{x,t)-n = 0, xedfl, (3.49) 

and the hyperbolic system (2.1), (2.2). The flux $ that we have to determine on the 
boundary dQ has the following form if we take into account (3.49): 
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Thus, we have to determine on the boundary dQ. For this purpose we use a comple- 
mentary numerical boundary condition as explained above: 

hiwK) ■ $ = hiwK) ■ J^n{wK), (3.51) 

where I3 is the left eigenvector corresponding to the positive eigenvalue A3 = m„ + c = c > 0. 
Solving equation (3.51) leads to the following value of the unknown component: 



2 



an ^ 2 ' 



K 



which determines completely the boundary flux (3.50). 



3.7.2. Generating boundary. Now let us consider a boundary where the total water depth 
is prescribed: 

H\q^^ = Ho{xs, t) > 0, Xs E d^l. 
Taking into account this information, the flux $ to be determined has the following form: 

fHoUr. 
Houun + ^H^n^ I . (3.52) 



Hence, we have to find u and v on the generating boundary dVt. The normal velocity will 
be immediately deduced from this information Un '■= un^ + vUy. 

Throughout this section we have assumed that the flow is "subsonic", i.e. \u ■ n\ < c. 
We could also consider the "supersonic" case, but physically this situation is rather exotic. 
Henceforth, we have one negative eigenvalue Ai = m„ — c, one positive A3 = m„ + c and 
A2 = Un can be in principle of any sign. Thus we have to consider two cases: m„ < and 
u„ > 0. In the first case we need a supplementary physical condition (on the tangential 
velocity to the boundary), in the second one we use a supplementary numerical condition: 

kiwx) ■ $ = kiwK) ■ J^n{wK)- 

Both lead to the same conclusion: Mrla^ = w,-!^, where Ur := uriy — mix is the tangential 
velocity. Computations similar to the previous section 3.7.1 lead to 



Substituting these expressions into (3.52) gives the boundary flux $. 



4. Numerical results 

Two kinds of numerical tests are presented. The first kind is a comparison with analytical 
solutions (or approximate analytical solutions): sections 4.1, 4.2, 4.4. This allows us to test 
the correctness and precision of the numerical scheme. The second kind is a comparison 
with results from laboratory experiments: section 4.3. This allows us to test the capacity of 
the code to reproduce actual events, and in particular to assess the validity of the nonlinear 
shallow-water equations for tsunami modeling. 
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Convergence rate in L norm 




h - average edge length 

Figure 13. Error of the numerical method in Loo norm. 

4.1. Convergence test. We begin the presentation of numerical tests by the simplest 
one - convergence test. We show the accuracy of the MUSCL scheme implementation. In 
order to do it, we solve numerically the following scalar linear advection equation 

— + Mo ■ Vw = 0, Mo e 
at 

with smooth^ initial conditions. Moreover, it has almost compact support in order to reduce 
the influence of boundary conditions. It is obvious that this equation will just translate 
the initial form in the direction uq. So, we have an analytical solution which can be used 
to quantify the numerical method error. On the other hand, to measure the convergence 
rate, we constructed a sequence of refined meshes. 

Figure 13 shows the error of the numerical method in L^o norm as a function of the mesh 
characteristic size. The slope of these curves represents an approximation to the theoretical 
convergence rate. On this plot, the blue curve corresponds to the first order upwind 
scheme while the other two (red and black) correspond to the MUSCL scheme with least- 
squares (see Section 3.6.4) and Green-Gauss (see Section 3.6.3) gradient reconstruction 
procedures respectively. One can see that the blue curve slope is equal approximatively 
to 0.97 which means first order convergence. The other two curves have almost the same 
slope equal to 1.90, indicating a second order convergence rate for the MUSCL scheme. In 
our implementation of the second-order scheme the least-squares reconstruction seems to 
give slightly more accurate results than the Green-Gauss procedure. 

The next figure represents the measured CPU time in seconds as a function of the mesh 
size. Obviously, this kind of data is extremely computer dependent but the qualitative 



We intentionally choose a smooth initial condition since the discontinuities can decrease the overall 
accuracy of the scheme. 
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Figure 14. CPU time for different finite volume schemes. 



behaviour is the same on all systems. On Figure 14 one can see that the "fastest" curve 
is the blue one (first order upwind scheme). Then we have two almost superimposed 
(black and red) curves referring to the second-order gradient reconstruction on variables. 
Here again one can notice that the least-squares method is slightly faster than the Green- 
Gauss procedure. On this figure we represented one more curve (the highest one) which 
corresponds to Green-Gauss gradient reconstruction on fluxes (it seems to be very natural 
in the context of the FVCF scheme explained in Section 3.1.1). Our numerical tests show 
that this method is more expensive from the computational point of view and we decided 
not to choose it. 

The next three sections deal with the validation of VOLNA against a set of benchmarks 
for tsunami modelling proposed at the Catalina 2004 workshop on long waves ; they are 
among the 6 benchmarks currently recommended by the United States National Oceanic 
and Atmospheric Administration (NOAA) for the evaluation of operational tsunami fore- 
casting models [SBT"'"07]. In order to help the reproductibility and comparison of numer- 
ical results, all the following test cases make use of publicly available data^. Although we 
present results for one mesh only for each benchmark, we have checked that the simulations 
converge as the mesh resolution increases. 

4.2. Tsunami run-up onto a plane beach. In this test case, we look at the runup 
of a tsunami wave over a plane sloping beach (of slope ^). The initial depression wave 
propagates leftwards (see figure 15). The result of the simulation is compared to an ana- 
lytical solution obtained by using the initial value problem technique of Carrier, Wu and 
Yeh [CWY03]. Note that the computational domain is approximately 50 km long, whereas 
the shoreline motion scale is 1 km ; hence, we choose to refine the mesh by a factor 10 near 

''' http : //nctr . pmel . noaa . gov/benchmark/ index . html 
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Figure 15. Catalina 1 benchmark — bathymetry and initial free surface, 
at the same vertical scale. 

the initial shoreline. The results presented here correspond to a resolution of 8 meters in 
the direction of propagation. Moreover, the initial free surface amplitude is a few meters 
high. Thus, due to the difference of spatial scales between the bathymetry, the initial free 
surface and the computational domain dimensions, the source term gUVh is very steep 
(see figure 15), which renders the use of a well-balanced scheme mandatory. This test 
case is one-dimensional. Since our code is two-dimensional, we implement it using a two 
dimensional computational domain, with translation invariance in the transverse direction. 

We can see on Figures 16 and 17 that the numerical results match pretty well the 
approximate analytical solution, especially near the shoreline location. This ensures the 
accuracy of the runup algorithm presented in section 3.3. 

4.3. Tsunami run-up onto a complex 3-dimensional beach. This experiment repro- 
duces at ^ scale the Monai valley tsunami, which struck the Island of Okushiri (Hokkaido, 
Japan) in 1993, in a 205 meters long wave tank. The computational domain reproduces the 
last 5 meters of the wave tank. The initial incident wave offshore is given by experimental 
data, and fed as a time dependent boundary condition. 

We compare the numerical results with the recorded data at three of the wave gages in- 
stalled in the wave tank : gages number 5, 7 and 9, of respective coordinates (4.521, 1.196), 
(4.521, 1.696) and (4.521, 2.196). This can be seen in Figures 19 - 21. The main wave (be- 
tween times 15 and 25 seconds) is very accurately described, at all gages we consider. 
Moreover, the maximal runup is adequately captured by the code. This value is extremely 
high, and occured in Monai Valley (the corresponding canyon in the experimental setup, 
along with the free surface at the moment of maximum elevation are shown in figure 22). 
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Figure 16. Catalina 1 benchmark — comparison between analytical (solid) 
and numerical (dashed) values for free surface at time 160 seconds. The gray 
line represents the beach, and the gray point the initial shoreline location. 

We obtain a maximal runup value of 8.05 centimeters, which corresponds to 32.2 meters 
put back at field scale. This is remarkably close to the measured value of 31.7 meters. 

Hence the numerical model is able to reproduce the laboratory scenario accurately, even 
without bottom friction modelling (and thus without any free parameter). In this realistic 
test case, the ability to refine the mesh near the zones of interest is a very nice asset. 

We performed another computation using the set-up described above. Namely, we solved 
equations (2.1), (2.2) completed by a conservation law for the wave energy. This system 
was recently proposed by Dutykh & Dias and we refer to [DD09a] for technical details 
and discussions. The total energy evolution is depicted on Figure 23. We represented 
two curves. The blue solid line corresponds to the solution of the augmented system of 
equations. The black dashed line refers to the total energy, estimated from conservative 
variables: 

1 1 
E^-pH\u\'' + -pgri', 

where p is the constant fluid density and rj = H — h is the free surface elevation with 
respect to the undisturbed water level. In complete accordance with results reported in 
[DD09a], the computed wave energy is not prone to numerical diffusion and has excellent 
monotonicity properties. Just at the end of the simulation one can notice a little decrease 
in both curves. In fact, it is induced by energy losses due to the wave run-up on the beach. 
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Figure 17. Catalina 1 benchmark — comparison between analytical (solid) 
and numerical (dashed) values for free surface at time 220 seconds. The gray 
line represents the beach, and the gray point the initial shoreline location. 

4.4. Tsunami generation and runup due to a 2-dimensional landslide. In this test 
case, a translating Gaussian shaped mass, initially at the shoreline, translates rightwards 




5 10 15 



Figure 18. Catalina 2 benchmark : initial free surface profile (left) 
bathymetry (right). 



time (s) 



36 



D. DUTYKH, R. PONCET, AND F. BIAS 




Figure 19. Catalina 2 benchmark — comparison between numerical 
and experimental data (red curve) at gage 5. 




Figure 20. Catalina 2 benchmark — comparison between numerical ijiffi^l|§) 
and experimental data (red curve) at gage 7. 



and creates a wave (see figure 24). The seafloor can be described by the following equation : 

h{x,t) = H{x) — ho{x,t), 



H{x) = xtan(/3), and ho{x,t) = 6exp 



3^/^ _ g 
5tan(/3) 



X being the direction of propagation. Here, S represent the maximum thickness of the 
sliding mass, /i the ratio between 6 and the horizontal length of the mass, and (3 the 
beach slope. Initially the submarine mass is partially underwater. Hence, this test case 
corresponds to a subaerial landslide. A sketch of this experiment can be seen in figure 24. 
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Figure 2 1 . Catalina 2 benchmark — comparison between numerical 
and experimental data (red curve) at gage 9. 




Figure 22. Maximal runup in Monai valley with a factor 3 magnification 
for the vertical scale. 

This benchmark is one- dimensional, but is implemented using a two-dimensional compu- 
tational domain, as in section 4.2. The numerical results we present are obtained as a one 
dimensional slice (which does not depend on the transverse variable). 

The result of the numerical simulation is compared to an analytical solution, computed as 
an approximate solution of the linear shallow water equations with a forcing term [LLS03]. 
In figure 25, we can see the comparison between the anal5d;ical and numerical wave surface 
profiles at times 16, 32 and 48 seconds, for P = 5.7°, S = 1 m, and /i = 0.01. These values 
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Figure 23. Computed and reconstructed wave energy profiles for Catalina 
2 bencfimark. 



ensure tliat we are in tlie linear regime of the shallow water equations (and thus that the 
analytical solution is a good approximation of the nonlinear equations solution). Hence, 
the comparison is meaningful. Good agreement is reached. 

We also performed wave energy computation for this test-case. To our knowledge, the 
energy evolution has not been shown yet for a landslide generated wave. Computation 
results are presented on Figures 26 and 27 for two cases /x = 0.01 and /i = 0.1. In the 
latter case the linear shallow water equations (LSWE) are not valid even on small time 
scales. See [LLS03] for more information. 




X (km) 



Figure 24. Sketch of the analytical landslide test case. The submarine 
mass displacement (in gray) and free surface (dashed line) have been mag- 
nified by a factor 1000. 
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Figure 25. Catalina 3 benchmark — free surface profiles at three differ- 
ent times — comparison between numerical results (dashed) and analytical 
formulas (solid lines), at times 16 (blue), 32 (red) and 48 (green) seconds. 
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Figure 26. Energy evolution with time for the 2-dimensional landslide test-case. 

Figure 26 shows the wave energy evolution with time. On Figure 27 we represented two 
trajectories in the energy phase-space {Ek,Ep), where is the kinetic energy and Ep is 
the potential one. We would like to point out the approximate energy repartition on the 
black curve (// = 0.1). 

4.5. Summary. Using different analytical benchmarks, we have validated all components 
of our code : accuracy and order of convergence of the numerical scheme (section 4.1), run 
up algorithm and treatment of steep depth fields (section 4.2), and time varying bathymetry 
in a conservative shallow- water framework (section 4.4). Moreover, we have shown the 
capability of the code to model realistic events, using an experimental benchmark (sec- 
tion 4.3). 
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Figure 27. Trajectories in the energy phase-space {Ek, Ep) for the 2D land- 
shde test-case. 

5. Conclusions and perspectives 

In the present article we provided a detailed description of the VOLNA code, designed 
for complete tsunami wave modelling. Namely, we are able to simulate the whole life-cycle 
of a tsunami from generation to inundation. Special attention was payed to the run-up 
algorithm described in Section 3.3. The overall performance test and validation were done 
in Section 4. 

The VOLNA code is operational and is able to run in complex and rapidly varying 
conditions. The use of unstructured meshes allows the handling of real coasts. Owing 
to the implementation of various types of boundary conditions, the code VOLNA can be 
coupled to other solvers and treat exclusively the zones where the NSWE are physically 
relevant. 

Some new results were presented concerning the energy of tsunami waves [DD09a]. In 
particular, we show the wave energy evolution for the Catalina 2 test case (run-up on a 
complex 3D beach) and a landslide generated tsunami (Catalina 3 benchmark problem). 

We are presently adding more physics to the VOLNA code: dissipative effects [DD07b, 
Dut09] (one could for example implement the dissipative terms from Bresch and Desjardins 
[BD03]) and dispersive effects. The current research activities focus on the development 
of robust finite volume solvers for dispersive wave models including the runup. The first 
attempt was made recently by D. Dutykh et al. (2011) [DKMll]. It was shown that 
dispersive effects might be beneficial for the description of the run-up of large, breaking 
waves. Finally a GPU acceleration of the VOLNA code is underway. Preliminary results 
can be found in [PVll]. 
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